#include<iostream>
#include<sstream>
#include<vector>
#include<set>
#include<math.h>
#include<gsl/gsl_randist.h>
#include "compute_alpha.hpp"


void readbins(){
  string line;
  int baf_count = 0;
  int rr_count = 0;
  while(getline(cin,line)){
    istringstream iss(line);
    feature_t feature;
    float baf;
    float lrr;
    iss>>feature.rs>>feature.chr>>feature.pos>>baf>>lrr;
    feature.maf = baf>.5?1-baf:baf;
    feature.rr = pow(2,lrr);
    if (feature.chr.compare("Y") && feature.chr.compare("XY") && feature.chr.compare("X") ){
      // compute the bafs
      for(int i=0;i<11;++i){
        sd_t & sd = sd_vec[i];
        if(feature.maf>=sd.mu_baf-CI_BAF && feature.maf<sd.mu_baf+CI_BAF){
          sd.rawdata_baf.push_back(feature.maf);
          ++baf_count;
          grand_sd_baf+=baf;
        }
        if(feature.rr>=sd.mu_rr-CI_RR && feature.maf<sd.mu_rr+CI_RR){
          sd.rawdata_rr.push_back(feature.rr);
          ++rr_count;
          grand_sd_rr+=feature.rr;
        }
      }
      if(feature.maf>=.5-CI_BAF && feature.maf<.5+CI_BAF){
        rr_offset+=(1.-feature.maf) * (feature.rr);
        het_weights+=(1.-feature.maf);
      }
      //cerr<<feature.chr<<" "<<feature.pos<<endl;
      bool newbin = false;
      if (bins.size()!=0){
        bin_t & bin = bins.at(bins.size()-1);
        unsigned int binsize = bin.size();
        if (binsize==0 || binsize<MAXBIN  && !bin.get_feature(binsize-1).chr.compare(feature.chr)  ){
          bin.add_feature(feature);
        }else{
          newbin = true;
        }
      }else{
        newbin = true;
      }
      if (newbin){
        bin_t bin;
        bin.add_feature(feature);
        bins.push_back(bin);
      }
    }
  }
  grand_sd_baf/=baf_count;
  grand_sd_rr/=rr_count;
  cerr<<"Grand BAF and RR SD: "<<grand_sd_baf<<" "<<grand_sd_rr<<endl;
  //rr_offset=pow(2,rr_offset/het_weights);
  rr_offset/=het_weights;
  cerr<<"RR offset is "<<rr_offset<<endl;
  //exit(0);
}
  
float normal_pdf(float x,float mu,float sigma){
   //float c = 2.506628;
  float p = 1/(sqrt(6.283185*sigma*sigma) * exp(.5 * pow((x-mu)/sigma,2)));
  return p;
}

void printbins(){
  float sd0_baf = sd_vec[10].sd_val_baf;
  float sd0_rr = sd_vec[10].sd_val_rr;
  for(int i=0;i<bins.size();++i){
    bin_t & bin = bins.at(i);
    float median = bin.get_median();
    // we will evaluate all 11 hypotheses and pick the max
    float maxLL = -1e20;
    float maxAlpha = -1;
    for(int k=0;k<11;++k){
      float alpha = sd_vec[k].alpha;
      //cerr<<"Working on alpha: "<<alpha<<endl;
      float LL = 0;
      for(int j=0;j<bin.size();++j){
        feature_t feature = bin.get_feature(j);
        float like_rr=normal_pdf(feature.rr+(1.-rr_offset),sd_vec[k].mu_rr,sd_vec[k].sd_val_rr);
        like_rr = 1;
        //cerr<<feature.rr<<" "<<rr_offset<<" "<<sd_vec[k].mu_rr<<" "<<sd_vec[k].sd_val_rr;
        float like_homo=normal_pdf(feature.maf,sd_vec[10].mu_baf,sd_vec[10].sd_val_baf);
        float like_het=normal_pdf(feature.maf,sd_vec[k].mu_baf,sd_vec[k].sd_val_baf);
        //if (like<1e-20) like = 1e-20;
     
        //cerr<<like_rr<<" "<<like_homo<<" "<<like_het<<endl;
        LL+=log(like_rr * ((1-HET_RATIO) * like_homo + HET_RATIO * like_het));
        //LL+=log((1-HET_RATIO)*normal_pdf(feature.maf,0,sd_vec[0].sd_val_baf)+(HET_RATIO)*normal_pdf(feature.maf,sd_vec[k].mu_baf,sd_vec[k].sd_val_baf));
      }
      cerr<<"Bin:"<<i<<" median:"<<median<<" LL for alpha "<<alpha<<": "<<LL<<endl;
      if (LL>maxLL){
        maxLL = LL;
        maxAlpha = alpha;
      }
    }
    cout<<i<<"\t"<<maxAlpha<<"\t"<<maxLL<<"\t"<<median<<endl;
    cerr<<i<<"\t"<<maxAlpha<<"\t"<<maxLL<<"\t"<<median<<endl;
    if (maxAlpha<0) exit(0);
  }
}

void init(){
  float alpha = 0;
  for(int i=0;i<11;++i){
    sd_t sd;
    sd.alpha = alpha;
    sd.mu_baf = (1.-alpha) * .5;
    sd.mu_rr = (1.-alpha) + alpha * (.32+1*.34); 
    cerr<<"Alpha: "<<alpha<<" murr "<<sd.mu_rr<<endl;
    sd.sd_val_baf = 0;
    sd.sd_val_rr = 0;
    alpha+=.1;
    sd_vec.push_back(sd);
  }
}

void compute_sd(){
  bool use_kai = true;
  if (use_kai){
    float sd0 = .016372;
    float sd_het_normal = .034982;
    for(int i=0;i<11;++i){
      //sd_vec[i].sd_val_baf = sd0+(sd_het_normal-sd0)*(.5-sd_vec[i].mu_baf)/.5;
      sd_vec[i].sd_val_baf = .015;
    }
    return;
  }
  bool use_grand = false;
  for(int i=0;i<11;++i){
    sd_t & sdobj = sd_vec.at(i);
    unsigned int n = sdobj.rawdata_baf.size();
    if(n){
      float mean = 0,sd=0;
      for(unsigned int j=0;j<n;++j){
        mean+=sdobj.rawdata_baf[j];
      }
      mean/=n;
      for(unsigned int j=0;j<n;++j){
        sd+=pow(sdobj.rawdata_baf[j]-mean,2);
      }
      sd = use_grand?grand_sd_baf:sqrt(sd/n);
      cerr<<"BAF STDDEV for alpha "<<sdobj.alpha<<": "<<sd<<endl;
      sd_vec[i].sd_val_baf = sd;
    }
    n = sdobj.rawdata_rr.size();
    if(n){
      float mean = 0,sd=0;
      for(unsigned int j=0;j<n;++j){
        mean+=sdobj.rawdata_rr[j];
      }
      mean/=n;
      for(unsigned int j=0;j<n;++j){
        sd+=pow(sdobj.rawdata_rr[j]-mean,2);
      }
      sd = use_grand?grand_sd_rr:sqrt(sd/n);
      cerr<<"RR STDDEV for alpha "<<sdobj.alpha<<": "<<sd<<endl;
      sd_vec[i].sd_val_rr = sd;
    }
  }
}

int main(int argc,char * argv[]){
  init();
  cerr<<"Initialized\n";
  readbins();  
  compute_sd();
  printbins();
  
}
